A dynamical equation for the distribution of a scalar advected by turbulence 
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A phenomenological model for the dissipation of scalar fluctuations due to the straining by the 
fluid motion is proposed in this letter. An explicit equation is obtained for the time evolution of 
the probability distribution function of a coarse-grained scalar concentration. The model relies on 
a self-convolution process. We first present this model in the Batchelor regime and then extend 
empirically our result to the turbulent case. The inclusion of this model in more general transport 
equations, including spatial gradients, is discussed in relation with 2D turbulence and stratified 
fiows. 
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The turbulent transport of tracers such as tempera- 
ture or salinity is an important issue for many applica- 
tions Pi . Available models usually provide a set of closed 
equations for the mean quantities and their second mo- 
ments (variance), and sometimes for the third and the 
fourth moments 2]. However, it may be useful to model 
the evolution of the whole probability distribution of the 
scalar. This problem has been discussed in the context of 
reactive flows 0- It could be also important to properly 
model the turbulent mixing of water masses in a stably 
stratified fluid, where sedimentation under gravity has 
an opposite effect for fluid particles heavier and lighter 
than the surrounding fluid. The case of vorticity in two- 
dimensional turbulence is also of interest. Indeed statis- 
tical mechanics of two-dimensional turbulence provides 
predictions for the final flow organisation depending on 
an initial distribution of vorticity values "^Q. This the- 
ory can be used as a starting point for non-equilibrium 
transport models 0, 0] , expressed in terms of the local 
probability density function of vorticity. Turbulent cas- 
cade however modifies this probability distribution in the 
presence of a small viscosity, with a dissipation of its fluc- 
tuations. In two-dimensional turbulence, this effect leads 
to a modification of the equilibrium state resulting from 
turbulent mixing 8] . We propose here a simple model for 
this cascade effect, which could be combined with trans- 
port equations in the presence of spatial gradients. 

Let us consider a scalar field cr(r,t) transported and 
conserved by the divergenceless turbulent motion of fiuid 
parcels. The effect of molecular diffusion is expected 
to smooth out the scalar field at the Batchelor's dif- 
fusive cut-off scale r^. We however prefer to consider 
the case of a purely advected scalar, with no diffu- 
sion, and introduce a local average at a given scale I, 
ai{v, t) — J Gi{r—r')a{r' ,t)dr' , obtained with a linear fil- 
tering operator Gi. This coarse-grained description cor- 
responds to a finite measurement resolution at a scale I 
larger than r^. The fine-grained probability distribution 
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of the scalar is preserved in time as the scalar value and 
the volume of each fluid parcel is conserved in the absence 
of diffusion. However the coarse-grained probability dis- 
tribution function (pdf) pi{a,t) is not preserved because 
fluctuations are transferred to scales smaller than the cut- 
off /. 

The problem is then to flnd a time evolution equation 
for this pdf. The result should of course depend on the 
properties of the turbulent field. In the usual Kolmogo- 
roff regime, the dissipation of scalar variance is equal to 
its cascade flux, independent of the cut-off scale. This 
flux is set by the energy and integral scale of the tur- 
bulence, which are described by transport equations in 
empirical turbulence models (like k-epsilon). 

Instead of a Kolmogoroff cascade, we shall here con- 
sider a random but smooth and persistent straining mo- 
tion, in which the velocity difference v(x + r,t) — v(x, t) 
is a linear function of the separation r |9| . This hypothe- 
sis holds whenever the kinetic energy spectrum is steeper 
than k~^. It could be an appropriate model for vorticity 
in two-dimensional turbulence. 

The straining is charaterized by the symmetric part 
of the strain tensor : — ^{djUi + diUj), while the 
antisymmetric part accounts for the rotation of a fluid 
element. As a consequence of the fluid incompressibility, 
the trace of S is zero and a basis exists on which this 
matrix is diagonal with two opposite eigenvalues ±s(i). 
The axes of an initial spherical blob corresponding to 
positive or negative eigenvalues will respectively grow or 
decrease. We assume that the angle between the posi- 
tive eigendirection of S and the iso-scalar lines evolves 
slowly compared to the straining rate s{t), i.e. the time 
for the eigendirections to rotate by 7r/2 is longer than 
the time for the scalar patterns to be strained from the 
integral scale L to the filtering scale I. With this hypoth- 
esis, our problem becomes locally one-dimensional: the 
fluid is composed of adjacent sheets of fluid uncorrelated 
with each other in one direction (which can be tilted with 
respect to the eigendirection of the strain matrix ) . 

After the time Ati/2 needed for the width of a strip 
to be divided by 2, the scalar field filtered at scale / be- 
comes the average of two realisations of the field at the 
previous time. The probabilities of scalar values in adja- 
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cent strips can be assumed independent, as they result of 
the straining of regions which were initially far apart, at 
a distance beyond the integral scale of the scalar. Thus 
the new probability distribution is the self-convolution of 
the previous one, describing the sum of the independent 
random variables, followed by a contraction by a factor 
2: p,(a,i + Ati/2) = 2 / pi{a' ,t)pi{<j - a' ,t)da' . 

The convolution can be transformed in a product of 
the Fourier transform of the pdf (characteristic function). 
If the scalar a has only positive values (which can be 
generally obtained by changing a in a ~ min{a)), it is 
more convenient to use the Laplace Transform pi{k) — 
J pi{a)e~'^"'da. The Laplace transform of the previous 
self-convolution relationship leads to pi{2n,t + Aii/2) = 
[pi{n,t)Y . Similarly, calling Atj_ the time to divide the 
thickness of a sheet of fluid by a factor n, the pdf at time 
t -I- Aij. will be an n*''-selfconvolution. In the spectral 
representation, this becomes the product of n identical 
characteristic functions [ij : 

p/(nK,t-f Aii/„) = [pi(K,i)]". (1) 

In order to get a differential equation in time, we now 
take the limit n = 1 + e with e = s{t)dt small in (Q. This 
yields: pi{k + eK,t + dt) = [pi{K,t)]^^'^ . Taking the limit 
e ^ 0, we can express pi{K + eK, t-\-dt) in terms of partial 
derivatives with respect to t and k: 



dtpi = s{t)[pi \npi - kOkPi] 



(2) 



One can check that the normalisation J pi(a,t)d(7 = 
pi{0,t) = 1 is preserved in time by this equation. The 
mean scalar value J crpi{a; t)da — — 9kP;(0, t) is also con- 
served (if it is initially defined). Note that if pi were 
standing for a Fourier instead of a Laplace transform, 
Pi{k) = J pi{a)e~'^'^'^ da , the same equation would be ob- 
tained for the norm of pi , with an additional equation for 
its phase cj): dt4>+ s{t)KdK<j) = s{t)(l>. 

The right hand term of equation Q describes the ef- 
fect of strain on the pdf of the scalar, which could be 
used with additional terms expressing scalar generation 
or spatial transport. Equation ([SJ itself can be analyt- 
ically solved (using the method of characteristics). For 
that purpose, we define the integral: 



/(i)=exp(/ s{t')dt') 
Jo 
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(such that s{t) ~ f'{t)/f{t)). We can easily show that 
a strip with initial width R{0) reaches a width R{t) — 
R{0)/f{t) at time t, so that f{t) is the reduction factor 
in the straining process. We can check that the result 
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is solution of jSJl. When f{t) is an integer, we recover 
the expression for the effect of n self-convolutions, in 
agreement with our initial assumption. 



The equation ^ and its solution Q describe the pro- 
cess of convergence to a Gaussian stated by the central 
limit theorem: at times goes on, scalar fluctuations ini- 
tially extending over more and more area become packed 
by the straining effect below the filtering scale. This can 
be be checked by the convergence to zero of all the cumu- 
lants beyond the second order one. The m*^ cumulant, 
defined as {cr"^)i.{t) = (— 9„)™ ln(pi(K, i))|^_Q, is readily 
obtained from Q), 
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The cumulant of order 2, equal to the variance, (cr^)^ = 

(cr^) — (cr)^, decays as l/f{t). This expresses the decay of 
the scalar variance by the cascade through the filtering 
scale. The relative value of the higher order cumulants is 

expressed as (o'™)^/(cr^}™^^ ^ f{t)^^™^'^^^\ so it decays 
in time, approaching a Gaussian, for which the cumu- 
lants with order larger than 2 are strictly equal to 0. 
Note however that if the first or second cumulants are 
not defined at time t = 0. the expression (^J converges 
to a Levy distribution another form of stable pdf, 
although the result |SJ for the cumulants would not be 
applicable. 

Equation |2Jl can be extented to the three dimensional 
case. The symmetric part of the deformation tensor may 
have 

• one negative eignevalue s{t): the problem remains 
one dimensionnal, but the fluid is now seen as a 
succession of adjacent iso-scalar planes 

• two negative eignevalues si{t), S2{t) : filaments are 
formed instead of sheets and the above description 
remains valid provided that s{t) — si{t) + S2{t). 

In the case of a scalar cascade in usual isotropic tur- 
bulence, our hypothesis of a uniform straining rate does 
not apply, but our approach can be still used in a more 
empirical way. We have seen indeed that the relative rate 
of decay of the scalar variance is equal to s{t), and this 
should be equal to the flux of scalar variance in the Kol- 
mogoroff cascade, independent of the cut-off scale. Then 
the model equation (jSJ can be applied to determine the 
evolution of the whole pdf. This approach ignores the 
fluctuations of s(t), which are known to generate inter- 
nal intermittence . Nevertheless, this can be a good 
model if cascade is of limited extent in wave numbers, 
or if larger sources of intermittence come from spatial 
gradients. 

The evolution of the pdf has been studied by from a 
numerical computation of isotropic turbulence, in which 
a passive scalar is introduced. The initial pdf was made of 
two symmetric peaks, and the convergence to a Gaussian 
pdf with decreasing variance was numerically observed. 
The corresponding result from our model is represented 
in figure (^b). The time scale is arbitrary, but s{t) could 
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istic function is 




a/<a> 



FIG. 1: Evolution of the shape of the scalar distribution 
as a sel-convolution process when the initial condition is a) 
strongly asymmetric (a decreasing exponential) b) symmetric 



be adjusted to fit the observed decay of the scalar vari- 
ance. Then the higher order cumulants can be calculated 
from the model. 

As a second example, we consider experiments per - 
formed by locally introducing a dye in a stirred flow |l4| . 
In that case the initial concentration is everywhere ex- 
cept in the dye streaks, so the corresponding pdf is very 
skewed. Later in time it progressively tends to a Gaus- 
sian as stirring proceeds. Villermaux and Duplat [l^ 
have provided a quantitative model of this evolution as 
an aggregation process of streaks of scalar, which leads 
to the following kinetic equation : 



(6) 



where f{t) is the integrated strain defined by relation Q 
like in our model. The first term accounts for the for- 
mation of scalar sheets, while the second one describes 
the decay of the concentration by the competing effects 
of strain and diffusion in scalar streaks. The solution 
of © approaches a sequence of gamma pdf at large 

time, 7(cr/(cr)) = r^T) ^^'^"^^^'^^ ' whose character- 
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The average concentration (cr) is a constant, while the 
exponent f(t) increases in time. The pdf is very skewed 
at the beginning and becomes more and more symmetric 
and narrow at time goes on. The solution Q has been 
found to be in good agreement with the experimental 
results obtained by introducing dye either in a steady 
stirring motion or in a turbulent pipe flow. 

Note that (7J) is in the form Q , so that it is also a solu- 
tion of our dynamical equation ^ . Therefore our model 
can also account for the experimental results. However 
in our case, it corresponds to a particular initial condi- 
tion, an exponential pdf (see figure I-a). Such an initial 
condition could result from non-homogeneous processes 
occuring near the dye injector. Other forms can be ob- 
tained from different initial conditions, although with the 
same qualitative behavior. By contrast, Villermaux and 
Duplat found that the family of gamma pdf is an attrac- 
tive solution 

111 

, so that it can be approached for a wider 
class of initial conditions. One can easily check that the 
dynamical equation © becomes identical to jSJ for long 
times, / oo. The change of the pdf in our model (j^J 
depends only on the straining rate s{t), while in (|SJ) it 
keeps track of the previous history through f{t), which 
may not be suitable in the presence of spatial fluxes. An- 
other difference with the model of Villermaux and Duplat 
[ll is that we consider concentration averaged on a small 
domain (in the absence of diffusion) instead of pointwise 
concentration. Distinguishing between the two models 
would require careful analysis of the experimental data. 

The effect of a turbulent diffusion in the presence of 
a mean scalar gradient has been studied by [T5| . who 
show that the pdf then develops exponential tails. This 
provides therefore a good rational for an initialisation 
of our cascade model by an exponential, leading then 
to solutions close to gamma pdf. In a steady regime 
sustained by a scalar gradient, the index / would then 
depends on the ratio of the cascade effect to the spatial 
fluxes, instead of time. 

Gamma pdf have been also proposed as a fit of density 
increments at meter scale in the ocean |l6j . Such distri- 
butions have been reproduced in numerical computations 
of internal wave breaking in a stably stratified fluid [l7i | . 
This was modeled as the effect of random steps in density 
due to mixing processes. However, such steps are not ob- 
served. Our approach could provide another justification. 
Indeed, a simple hypothesis of statistical equilibrium in 
the gravity field leads to an exponential pdf depending on 
density a and vertical height z. This would provide the 
Gamma pdf through the self-convolution process. How- 
ever the link between the statistics of this coarse-grained 
concentration at scale I and those of vertical density in- 
crements at scale / would need to be claridied. 

To conclude, the main result of this letter is the con- 
tinuous process Q that accounts for the temporal evolu- 
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tion of the coarse-grained scalar pdf : the probe of width 
I "sees" structures coming from larger and larger scales, 
which implies the self-convolution of the pdf. . The ef- 
ficiency of the fluid motion to drive the scalar pdf to a 
form given by the central limit theorem depends of the 
rate of strain s{t). It does not take into account inter- 



mittency by assuming that s{t) is homogeneous in space. 
Spatial gradients could in some cases be the main raison 
of non-gaussianity of the scalar distribution. In that re- 
spect, a fuller study using this result in the context of 
statistical mechanics of the mixing of a stably stratified 
fluid is in preparation and will be reported elsewhere. 
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